#ifndef UnstructMG_POINT_GS_HPP
#define UnstructMG_POINT_GS_HPP

#include "par_CSRMatrix.hpp"

// data_t是数据存储的精度，calc_t是算子作用时的精度（对应向量的精度）
template<typename idx_t, typename data_t, typename setup_t, typename calc_t, int dof>
class PointGS : public Solver<idx_t, setup_t, calc_t, dof> {
public:
    typedef enum {BACK_FORW, BACKWARD, SYMMETRIC, FORWARD, FORW_BACK, RED_BLACK} SCAN_TYPE;
    // 对称GS：0 for sym, 1 for forward, -1 backward
    SCAN_TYPE scan_type = SYMMETRIC;

    PointGS();
    PointGS(SCAN_TYPE type);
    ~PointGS();
    void SetOperator(const Operator<idx_t, setup_t, dof> & op);
    void SetScanType(SCAN_TYPE type);

protected:
    mutable bool last_time_forward = false;

    CSRMatrix<idx_t, data_t, calc_t, dof> sep_L;// separate L part of diagonal matrix
    CSRMatrix<idx_t, data_t, calc_t, dof> sep_U;// separate U part of diagonal matrix
    std::vector<calc_t>                   inv_D;// separate D part of diagonal matrix
    
    const par_CSRMatrix<idx_t, setup_t, setup_t, dof> * A_calc = nullptr;// 仅用于提供通信句柄
    bool owns_offd = false;
    const CSRMatrix<idx_t, data_t, calc_t, dof> * pgs_offd = nullptr;
    calc_t * offd_x = nullptr;
    calc_t * send_buf = nullptr;
    calc_t * tmp_rhs = nullptr;

    idx_t * map_old2new = nullptr;
    idx_t * map_new2old = nullptr;

    void separate_LU(const CSRMatrix<idx_t, setup_t, setup_t, dof> & A);    
    
    void truncate();
    // 近似求解一个（残差）方程，以b为右端向量，返回x为近似解
    void Mult_ipl(const par_Vector<idx_t, calc_t, dof> & b, par_Vector<idx_t, calc_t, dof> & x) const;
    
private:
    void wrapper_forward (const calc_t * rhs, calc_t * sol) const;
    void wrapper_backward(const calc_t * rhs, calc_t * sol) const; 

    void forward (const calc_t * rhs, calc_t * sol) const;
    void backward(const calc_t * rhs, calc_t * sol) const;
};

typedef PointGS<int, double, double, double, 1> PointGS_I32All64_DOF1;
typedef PointGS<int, double, double, double, 2> PointGS_I32All64_DOF2;
typedef PointGS<int, double, double, double, 3> PointGS_I32All64_DOF3;
typedef PointGS<int, double, double, double, 4> PointGS_I32All64_DOF4;

typedef PointGS<int, float , double, double, 1> PointGS_I32K64P64D32_DOF1;
typedef PointGS<int, float , double, double, 2> PointGS_I32K64P64D32_DOF2;
typedef PointGS<int, float , double, double, 3> PointGS_I32K64P64D32_DOF3;
typedef PointGS<int, float , double, double, 4> PointGS_I32K64P64D32_DOF4;

typedef PointGS<int, float , double, float , 1> PointGS_I32K64P32D32_DOF1;
typedef PointGS<int, float , double, float , 2> PointGS_I32K64P32D32_DOF2;
typedef PointGS<int, float , double, float , 3> PointGS_I32K64P32D32_DOF3;
typedef PointGS<int, float , double, float , 4> PointGS_I32K64P32D32_DOF4;

typedef PointGS<int, float , float , float , 1> PointGS_I32All32_DOF1;
typedef PointGS<int, float , float , float , 2> PointGS_I32All32_DOF2;
typedef PointGS<int, float , float , float , 3> PointGS_I32All32_DOF3;
typedef PointGS<int, float , float , float , 4> PointGS_I32All32_DOF4;
#ifdef USE_FP16
typedef PointGS<int, __fp16, double, float , 1> PointGS_I32K64P32D16_DOF1;
typedef PointGS<int, __fp16, double, float , 2> PointGS_I32K64P32D16_DOF2;
typedef PointGS<int, __fp16, double, float , 3> PointGS_I32K64P32D16_DOF3;
typedef PointGS<int, __fp16, double, float , 4> PointGS_I32K64P32D16_DOF4;

typedef PointGS<int, __fp16, float , float , 1> PointGS_I32K32P32D16_DOF1;
typedef PointGS<int, __fp16, float , float , 2> PointGS_I32K32P32D16_DOF2;
typedef PointGS<int, __fp16, float , float , 3> PointGS_I32K32P32D16_DOF3;
typedef PointGS<int, __fp16, float , float , 4> PointGS_I32K32P32D16_DOF4;
#endif
#endif